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Abstract 

Neural systems encode external stimulation information in the firing patterns of individual neurons. This paper deciphers how 
neural systems represent and transmit acupuncture information. First, the firings of the individual neuron are transformed into 
point process spike trains. Then a Bayesian approach to this problem based on state-space representations of point processes is 
described. We use this method to decode different types of manual acupuncture (MA), i.e. the stimulus oscillograms of the MAs 
are reconstructed from neural spike trains. These results are proved to be significant for neural data processing based on 
statistical analysis framework. 
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Introduction 

The relationship between neural spiking activity and a stimulus or behavior is the focus of neuro physiological 
investigation^]. Describing the way a neuron represents information about a stimulus or behavior is the encoding 
problem. The dual problem of reproducing the stimulus or behavior from neural spike trains is the decoding 
problem. The ability to reconstruct a signal from observed neural activity provides a way of verifying that explicit 
information about the external world is represented in particular neural spike trains [2-3]. 

Since the spikes from a single neuron have stereotyped waveforms, neural representations of external stimulus 
signals are based on the frequency and timing of spike events. In addition, the neural spike trains generated in 
response to identical stimuli applied on separate occasions will not be identical and exhibit substantial variation 
even they may share common statistical features. These two fundamental facts suggest that neural spike train data 
might be effectively analyzed with the theory of point processes [4-6]. In this case, the temporal point processes are 
used to represent the times of the spiking events and to express the probability distribution of specific sequences of 
spike times. 

In this article, we review the state-space approach for solving the neural decoding problem. A state-space 
estimation and inference framework is constructed by writing state models that describe the evolution of the 
stochastic signals to estimate, and intensity models for neurons that define the probability density of observing a 
particular sequence of spike times. 

Here we establish the relationship between acupuncture stimulus oscillogram[7] and electrical signals of neuron 
through state-space model and reconstruct acupuncture signals of different types. This is a preliminary exploration 
that why MA with different types have different effects. 

Experimental Materials and Methods 

Description of Acupuncture Manipulation 

The lifting-thrusting MA manipulation method contains three different types of manipulations based on the 
duration of lifting phase and thrusting phase, and they are lifting-thrusting reinforcing (' I B'), lifting-thrusting 
reducing ('TX') and uniform lifting-thrusting ('TC') manipulations respectively. When we carry out the MA 
manipulation, the depth of the acupuncture needle in the skin is real- timely recorded as an oscillogram [7], see Fig. 1. 
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FIG. 1 SCHEMATIC DIAGRAM ABOUT THE OSCILLOGRAMS OF 'TB', 'TX' AND 'TC' MA 
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FIG. 2 SPIKE-SORTING OF SIGNALS 


Acupuncture Neural Data Preprocessing 

The neural spike activity evoked by acupuncture was recorded as continuous wide-band data, so spike detection 
and sorting should be accomplished firstly [7] to obtain spike train of individual neuron. After band-filtering the 
data between 300 and 3000 Hz, the spikes were detected by setting the threshold at 

\ \x\ i 

Thr = 3cr, cr = median { — — — > (1) 

[ 0.6745 J 

where x is signal after band-filtering and a is the estimate of the standard deviation of the noise. The time series 
after band-filtering are shown in Fig. 2. 

After the detection of spikes, the features of the spike shapes was extracted by the wavelet transform as inputs to 
the clustering algorithm. Each acupuncture stimulus evokes an obvious bursting, which contains the main 
information of acupuncture encoded by neuron. Then the most obvious class of spike train is transformed to point 
process. 

Neural Spiking Model and Neural Decoding 

Conditional Intensity Function 

A point process is defined as a set of discrete events that occur in continuous time. For a neural spike train this 
would be the set of individual spike times[8]. Assume that the spiking activity of a single neuron is recorded for an 
observation interval (0, T ] during a neurophysiological experiment. Let 0<u l <---<u J <---<u J <T be the set of J 

spike times from the neuron, i.e. the point process observations. For / e ( 0, 7’] , we define N(t) as the number of 
spikes counted in the time interval (0, t\ . Then the sample path of the spike times from the neuron is defined as the 
event N 0t = jo < u x < u 2 ,...,Uj < t f|7V(f) = / j with j < J and it jumps 1 at the spike times and is constant otherwise. 
The sample path tracks the location and number of spikes in (0, t\ and therefore contains all the information of 
spike times[9]. 

A stochastic neural point process can be completely characterized by its conditional intensity function[8], which is 
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defined as 

^|// f ) = lim PW + A) -^> = 1| ^ ) (2) 

A->0 A 

where H t is the neuron's spiking history up to and includs time t and other relevant covariates. The conditional 
intensity function is a history-dependent rate function that generalizes the definition of the rate function of a 
Poisson process[8]. 

Using the conditional intensity function A(t | H t ) , the joint probability distribution of a neural spike train {w ; }" , in 
an interval (0,r] is expressed as 


pmi i) = 


n Mu, 

i=l 


H(u,)) 


exp 


I H(u))) 


(3) 


It is obviously that the first term on the right-hand side of (3) can be viewed as the probability density associated 
with observing action potentials at the actual spike times and the second term as the probability of not observing 
any other spikes anywhere in the interval. 


Neural Point Process Model 


Multiple factors simultaneously affect spiking activity of single neurons. The class of conditional intensity models 
that we use allows the analysis of the simultaneous effects of extrinsic covariates, internal spiking history, and 
concurrent neural ensemble activity. That is the conditional intensity processes with which we work can be written 
generally as 


T(t|//,) = /(hx [0 , (] ,/V [0;() ,{/V [ ‘ 0?) } c c =1 ) 


(4) 


where x [0;] represents the values of a set of external covariates up to and including the present time, N, 0 t 
represents the neuron's spiking history, and VUj jy , represents the firing histories of other C neurons. 


Here, we model a simple form of the conditional intensity function [9] in terms of the state process as 

A(k A) = exp (p + J3x k ) 


(5) 


where p is the log of the background firing rate and // is its gain parameter that governs how much the latent 
process modulates the firing rate of this neuron. 


Assessing Model Goodness- of -fit 

Before making an inference about a particular neural system from a statistical model, it is crucial to measure 
agreement between the model and the spike train data series, that is, to evaluate goodness-of-fit [9]. The time- 
rescaling theorem is a well-known result in probability theory to assess goodness-of-fit, and provides a direct 
means of measuring agreement between point process and a probability model intended to describe its stochastic 
structure [10]. The Kolmogorov-Smirnov test rejects the null-hypothetical model if any of the plotted points lie 
outside these bands. 


Neural Decoding Analysis by State Estimation 

Based on the discrete-time framework, the Bayesian decoding process can be constructed by using nonlinear 
recursive filter and fixed interval smoothing algorithm[9]. 


Decoding of Acupuncture Signals 

State-Space Model of Acupuncture Electrical Signals 

Now we illustrate the state-space framework for acupuncture electrical signals and solve the neural decoding 
problem. Neurons in spinal dorsal root ganglion fire preferentially when it response to the acupuncture stimulus. 
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The observation suggests that the acupuncture oscillogram can be reconstructed from neural spiking activity. We 
take the conditional intensity function of the neuron to be 

^(**) = exp (/i + p-x k ) (22) 

where p a □ 1 sets the baseline firing rate, and /? e □ 1 is the response strength of the neuron to acupuncture. 
Because the neurons are responsible for the generation and conduction of the electrical signal evoked by 
acupuncture in the initial stages, we set up this model under the assumption that the neuron has firing properties 
that do not depend on its own spike history or that of any other neuron. Therefore this activity comprises an 
inhomogeneous Poisson process. The observation interval is partitioned into uniform bins with size of A . The state 
evolution equation is taken to be 

X k = px, t_i + s k (23) 

where p is a correlation coefficient, s k is a Gaussian random variable with mean zero and variance a] . 
Reconstructing Acupuncture Oscillogram 

We applied the Bayesian decoding algorithm to data recorded from the spinal dorsal root ganglion of experiment 
rats, on which different types of MA ('TB','TX' and 'TC') are taken. The true acupuncture oscillograms of different 
types ('TB','TX' and 'TC' MA) are also simulated based on the acupuncture stimulus times. The whole time interval 
for each record is T = 20 seconds. The resolution of the spike trains is 1 msec, and the state equations are updated 
at 1 msec. 

Similar to the analysis of signals generated by simulation, acupuncture oscillograms are also reconstructed from 
the neural spike trains by Bayesian decoding algorithm. First, we estimate the parameters in equations by 
maximum likelihood estimation. The parameter estimates from the EM algorithm are shown in Table 1. 


TABLE 1 THE PARAMETER ESTIMATES FOR ACUPUNCTURE SIGNALS. 
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FIG. 3 THE ESTIMATED OSCILLOGRAMS 
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4 THE K-S PLOTS OF THE MODEL ESTIMATE 
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After the parameters were estimated, acupuncture oscillograms were reconstructed from neural spike trains. Fig. 6 
shows a fragment of the estimated states of three different types of MA respectively. For the 'TB' MA 
manipulations, the oscillogram has the bigger duration of ascending branch (tl) and smaller duration of 
descending branch (t2), as seen from Fig. 6(a); for the 'TX' MA manipulations, the oscillogram has the smaller 
duration of ascending branch (tl) and bigger duration of descending branch (t2), as seen from Fig. 6(b); the 
durations of ascending branch (tl) and descending branch (t2) are approximately the same in the oscillogram of 
TC' MA manipulation, as seen from Fig. 6(c). Overall, decoded acupuncture oscillograms captured mostly 
fluctuations as the simulated state shown and showed the basic characteristics of each MA manipulations. The K-S 
plots are shown in Fig. 7 to assess model goodness-of-fit. For 'TB' MA, the models lie within the 95% confidence 
limits, indicating close agreement between the overall model fit and the data. For both 'TX' and 'TC' MA, the fit of 
model lies within the confidence limits except at the small quantiles. 

Conclusions 

We have introduced Bayesian decoding algorithm to the studies of acupuncture neural electrical signals, i.e. we 
construct the relationship between different types of MA and the corresponding neural spike trains. The neural 
spike trains can be described by point process modeling based on the conditional intensity function. To reconstruct 
acupuncture oscillogram from point process spike trains we have used state-space models which are well suited 
for describing these changes in neural firing properties. The present study shows that the encoding of different 
acupuncture manipulations have their own features. Different MA oscillogram can be reconstructed from the 
neural spike trains by decoding algorithm even in a single trial. These results imply that different encoding 
features are related to different effects which different acupuncture manipulations have. 
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